#include <iostream>
	using std::cerr;
#include <limits>
	using std::numeric_limits;
#include <iterator>
	using std::back_inserter;
#include <string>
	using std::string;
#include <algorithm>
	using std::max;
	using std::min;
#include <vector>
	using std::vector;


#include <boost/range/algorithm.hpp>
	using boost::range::lower_bound;
#include <boost/range/adaptor/map.hpp>
	using boost::adaptors::map_keys;
#include <boost/range/algorithm/copy.hpp>
	using boost::range::copy;
#include <boost/filesystem.hpp>
	namespace fs = boost::filesystem;
#include <boost/filesystem/fstream.hpp>
	//using fs::ifstream

#include <boost/spirit/include/karma.hpp>

#include <boost/format.hpp>
	using boost::format;
#include <boost/tr1/tuple.hpp>
//#include <boost/assign.hpp>
//#include <boost/assign/std/vector.hpp> 
#include <boost/assign/list_of.hpp> 
	using namespace boost::assign; // bring 'list_of()' into scope

#include "discrete_intervals_per_contig.h"
#include "parse_genomic_intervals.h"
#include "sort_by_embedded_numbers.h"


// N.B. uses open-close coordinates

namespace discrete_intervals_per_contig
{


//
//	special magic intervals for out of bounds
//
t_intervals out_of_bounds_intervals = list_of(t_interval(0,0))(t_interval(numeric_limits<unsigned>::max(),numeric_limits<unsigned>::max()));

//________________________________________________________________________________________
//
//  c'tor
// 
//________________________________________________________________________________________
t_discrete_intervals::t_discrete_intervals()	:
		i_prev_interval(out_of_bounds_intervals.begin()),
		i_next_interval(out_of_bounds_intervals.begin()),
		curr_contig_intervals(intervals_by_contig.end())
{}

void t_discrete_intervals::init()
{
	//	
	//  Sort intervals
	// 
	using boost::range::sort;
	for (i_intervals_by_contig ii = intervals_by_contig.begin(), end = intervals_by_contig.end(); ii != end; ++ii)
		sort((*ii).second, t_sort_by_beg());
}


//________________________________________________________________________________________
//
//  set_curr_contig
// 
//________________________________________________________________________________________
bool t_discrete_intervals::set_curr_contig(const string& test_contig)
{
	// 	
	// Returns false if contig is absent
	// 
	i_intervals_by_contig contig_intervals = intervals_by_contig.find(test_contig);
	

	// do nothing if missing contig
	if (contig_intervals == intervals_by_contig.end())
		return false;

	if (contig_intervals->second.size() == 0)
		return false;

	// 
	//  Put everything in a new sane state
	// 
	curr_contig = test_contig;
	curr_contig_intervals = contig_intervals;
	i_prev_interval = curr_contig_intervals->second.begin();
	i_next_interval = i_prev_interval;
	return true;

}


//________________________________________________________________________________________
//
//  contains
// 
//________________________________________________________________________________________
bool t_discrete_intervals::contains(string test_contig, unsigned test_pos)
{
	#ifdef DEBUG
	cerr << "    Check if contains " << test_pos << "?\n";
	#endif
	//
	//	new contig 
	// 
	if (curr_contig != test_contig)
	{	
		#ifdef DEBUG
		cerr << "    Different Contig from previous test\n";
		#endif
		// 
		//  Missing contig
		if (!set_curr_contig(test_contig))
			return false;
		//
		//	Get next interval 
		// 
		get_flanking_intervals(test_pos);
		return (i_next_interval == i_prev_interval);
	}

	// within flanking intervals. Just test if within a single interval
	if (within_flanking_interval_pair(test_pos))
	{
		#ifdef DEBUG
		cerr << format("    Same interval as previous test: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
		#endif
		return (i_next_interval == i_prev_interval);
	}
	else
	{
		#ifdef DEBUG
		cerr << format("    Different Interval: Previous flanking intervals = (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
		#endif

	}


	// 
	// increment to next interval
	//  	Takes two increments if already within the bounds of an interval
	//  	Otherwise, if between intervals, requires a single interval
	if (i_next_interval == i_prev_interval)
	{
		increment_flanking_interval_pair();
		if (within_flanking_interval_pair(test_pos))
		{
			#ifdef DEBUG
			cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
			#endif
			return (i_next_interval == i_prev_interval);
		}
		increment_flanking_interval_pair();
		if (within_flanking_interval_pair(test_pos))
		{
			#ifdef DEBUG
			cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
			#endif
			return (i_next_interval == i_prev_interval);
		}
	}
	else
	{
		increment_flanking_interval_pair();
		if (within_flanking_interval_pair(test_pos))
		{
			#ifdef DEBUG
			cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
			#endif
			return (i_next_interval == i_prev_interval);
		}
	}

	get_flanking_intervals(test_pos);
	#ifdef DEBUG
	cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
	#endif
	return (i_next_interval == i_prev_interval);
}


//________________________________________________________________________________________
//
//  get_flanking_intervals
// 
//________________________________________________________________________________________
void t_discrete_intervals::get_flanking_intervals(unsigned test_pos)
{
	t_intervals& intervals = curr_contig_intervals->second;
	 
	// 
	//	check not past the end 
	// 
	i_next_interval = lower_bound(intervals, t_interval(test_pos, test_pos));


	// If test_pos is before first interval: 
	// 
	//  	prev_interval => out_of_bounds_intervals[0]
	//  	next_interval => intervals[0]
	if (i_next_interval == intervals.begin() && test_pos < i_next_interval->beg)
	{
		#ifdef DEBUG
		cerr << "    get_flanking_intervals: Before any intervals\n";
		#endif
		i_prev_interval = out_of_bounds_intervals.begin();
		return;
	}


	// after last interval
	// 
	//  	next_interval => out_of_bounds_intervals[1]
	//  	prev_interval => intervals[-1]
	// 
	if (i_next_interval == intervals.end())
	{
		#ifdef DEBUG
		cerr << "    get_flanking_intervals: Beyond any intervals\n";
		#endif
		i_prev_interval = intervals.begin() + intervals.size() - 1;
		i_next_interval = out_of_bounds_intervals.begin() + 1;
		return;
	}

	// entirely contained:
	// both flanking intervals point to the same thing
	// 
	if (test_pos >= i_next_interval->beg) // && test_pos < i_next_interval->end)
	{
		#ifdef DEBUG
		cerr << "    get_flanking_intervals: Contained within a single interval\n";
		#endif
		i_prev_interval = i_next_interval;
		return ;
	}

	// 
	// before next
	// 
	i_prev_interval = i_next_interval - 1;
    #ifdef DEBUG
	cerr << "    get_flanking_intervals: Between two intervals\n";
    #endif
	return;
}

//________________________________________________________________________________________
//
//  any_closest
// 
//________________________________________________________________________________________
void t_discrete_intervals::increment_flanking_interval_pair()
{

	// can't increment any more
	if (i_next_interval == out_of_bounds_intervals.begin() + 1)
	{
		#ifdef DEBUG
		cerr << "    increment_flanking_interval_pair: Already beyond any intervals\n";
		#endif
		return;
	}

	// if both point to the same interval, increment the next iterator so that they once
	//	again flank
	if (i_prev_interval == i_next_interval)
	{
		++i_next_interval;
		t_intervals& intervals = curr_contig_intervals->second;
		if (i_next_interval == intervals.end())
		{
			#ifdef DEBUG
			cerr << "    increment_flanking_interval_pair: Beyond any intervals after increment\n";
			#endif
			i_next_interval = out_of_bounds_intervals.begin() + 1;
		}
		else
		{
			#ifdef DEBUG
			cerr << "    increment_flanking_interval_pair: Between intervals after increment\n";
			#endif
		}
		return;
	}

	//  Have both point within the next interval
	i_prev_interval = i_next_interval;
    #ifdef DEBUG
	cerr << "    increment_flanking_interval_pair: Within single interval after increment\n";
    #endif
	return;

}
//________________________________________________________________________________________
//
//  within_flanking_interval_pair
// 
//________________________________________________________________________________________
bool t_discrete_intervals::within_flanking_interval_pair(unsigned test_pos, unsigned& dist, i_intervals& closest_interval) const
{
	if (i_prev_interval == i_next_interval)
	{
		// 
		// Entire contained within single interval
		// 
		// 		Both interval iterators point to the same interval
		// 
		if (test_pos >= i_prev_interval->beg && 
			test_pos <  i_prev_interval->end  )
		{
			dist = 0;
			closest_interval = i_prev_interval;
			return true;
		}
		return false;
	}
	else
	{
		if (test_pos >= i_prev_interval->end && 
			test_pos <  i_next_interval->beg  )
		{
			unsigned dist1 = test_pos - i_prev_interval->end + 1;
			unsigned dist2 = i_next_interval->beg - test_pos;
			
			//
			//  Choose next interval if
			//  	1) NOT out_of_bounds end
			//  	2) shorter dist / out_of_bounds beg
			// 
			if ((i_next_interval != out_of_bounds_intervals.begin() + 1) &&
				(dist1 > dist2										||
				 i_prev_interval == out_of_bounds_intervals.begin() 
				 ))
			{
				dist = dist2;
				closest_interval = i_next_interval;
			}
			else
			{
				dist = dist1;
				closest_interval = i_prev_interval;
			}
			return true;
		}
	}
	return false;
}

bool t_discrete_intervals::within_flanking_interval_pair(unsigned test_pos) const
{
	if (i_prev_interval == i_next_interval)
	{
		// 
		// Entire contained within single interval
		// 
		// 		Both interval iterators point to the same interval
		// 
		if (test_pos >= i_prev_interval->beg && 
			test_pos <  i_prev_interval->end  )
			return true;
		return false;
	}
	else
	{
		if (test_pos >= i_prev_interval->end && 
			test_pos <  i_next_interval->beg  )
			return true;
	}
	return false;
}


//________________________________________________________________________________________
//
//  any_closest
// 
//  	only return false if contig is absent
// 
//________________________________________________________________________________________
bool t_discrete_intervals::any_closest(string test_contig, unsigned test_pos, unsigned& dist, i_intervals& closest_interval)
{
	//
	//	new contig 
	// 
	if (curr_contig != test_contig)
	{	
		// 
		//  Missing contig
		if (!set_curr_contig(test_contig))
			return false;
		//
		//	Get next interval 
		// 
		get_flanking_intervals(test_pos);
		within_flanking_interval_pair(test_pos, dist, closest_interval);
		return true;
	}

	// try current flanking interval pair
	if (within_flanking_interval_pair(test_pos, dist, closest_interval))
	{
		#ifdef DEBUG
		cerr << format("    Same interval as previous test: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
		#endif
		return true;
	}
	#ifdef DEBUG
	cerr << format("    Different Interval: Previous flanking intervals = (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
	#endif

	// or next one
	increment_flanking_interval_pair();
	if (within_flanking_interval_pair(test_pos, dist, closest_interval))
	{
		#ifdef DEBUG
		cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
		#endif
		return true;
	}

	// or next one
	increment_flanking_interval_pair();
	if (within_flanking_interval_pair(test_pos, dist, closest_interval))
	{
		#ifdef DEBUG
		cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
		#endif
		return true;
	}

	//
	//	OK. do proper search
	// 
	get_flanking_intervals(test_pos);
	within_flanking_interval_pair(test_pos, dist, closest_interval);
	#ifdef DEBUG
	cerr << format("    New Interval: (%1$d-%2$d) - (%3$d-%4$d)\n") % i_prev_interval->beg  % i_prev_interval->end  % i_next_interval->beg  % i_next_interval->end;
	#endif
	return true;
}

//________________________________________________________________________________________
//
//		expand_by
// 
//________________________________________________________________________________________
void t_discrete_intervals::expand_by(unsigned by)
{
	//
	//  iterate through each set of intevals
	// 
	for (i_intervals_by_contig ii = intervals_by_contig.begin(), ee = intervals_by_contig.end(); ii != ee; ++ii)
	{
		// current intervals
		t_intervals& intervals = ii->second;
		if (intervals.size() < 1)
			continue;

		i_intervals src = intervals.begin();
		i_intervals end = intervals.end();

		#ifdef DEBUG
		std::cerr << format("Expand Orig : (%1$d-%2$d) by %3$d [%4$s]\n") % src->beg % src->end % by % ii->first;
		#endif
		while (src != end)
		{
			// make sure no overflow
			src->beg = max(by, src->beg) - by;
			src->end = min(numeric_limits<unsigned>::max() - by, src->end) + by;
			++src;
		}
		#ifdef DEBUG
		std::cerr << format("New Interval: (%1$d-%2$d)\n") % intervals.begin()->beg % intervals.begin()->end;
		#endif
	}
	combine_overlapping();

}

//________________________________________________________________________________________
//
//  combine_close_neighbours
// 
//  	Merge neighbours that are less than a certain   
//  	distance apart.
// 
//  	side effect: sort intervals
//________________________________________________________________________________________
void t_discrete_intervals::combine_close_neighbours(unsigned min_distance)
{
	// 
	// invalidate cached intervals
	// 
	i_prev_interval			=out_of_bounds_intervals.begin();
	i_next_interval			=out_of_bounds_intervals.begin();
	curr_contig_intervals	=intervals_by_contig.end();
	curr_contig				= string();

	//
	//  iterate through each set of intevals
	// 
	for (i_intervals_by_contig ii = intervals_by_contig.begin(), ee = intervals_by_contig.end(); ii != ee; ++ii)
	{
		// current intervals
		t_intervals& intervals = ii->second;
		if (intervals.size() < 2)
			continue;

		//
		//  sort
		// 
		using boost::range::sort;
		sort(intervals, t_sort_by_beg());

		i_intervals dst = intervals.begin();
		i_intervals src = intervals.begin() + 1;
		i_intervals end = intervals.end();

		//
		//  merge in place, keeping id of first interval in set of overlapping
		// 
		while (src != end)
		{
			if (src->beg  < min_distance + dst->end + 1)
			{
				// merge into current dst
				dst->end = max(src->end, dst->end);
			}
			else
			{
				// keep current dst, inc and copy src->dst
				++dst;
				*dst = *src;
			}
			++src;
		}

		++ dst;
		intervals.erase(dst, end);
	}
}


//________________________________________________________________________________________
//
//  read_intervals
// 
//________________________________________________________________________________________
void read_intervals(std::istream& f_intervals, t_discrete_intervals& intervals, bool inclusive_coordinates)
{
	t_intervals_by_contig& intervals_by_contig = intervals.intervals_by_contig;

	// adjust values to convert from inclusive to open close coordinates
	unsigned adjust_value = inclusive_coordinates ? 1 : 0;
	unsigned cnt_line = 0;

	sp_genomic_interval_grammar genomic_interval_grammar = create_genomic_interval_grammar(cnt_line);

	string str;
	// ignore first line
	while (std::getline(f_intervals, str))
	{
		cnt_line ++;
		if (!str.length() || str[0] == '#')
			continue;

		t_genomic_interval interval;
		if (!parse_genomic_interval(genomic_interval_grammar, str, interval))
			break; 

		// convert to open / closed by adding 1
		intervals_by_contig[interval.contig].push_back(t_interval(interval.beg - adjust_value, interval.end));
	}
	intervals.init();
}

//________________________________________________________________________________________
//
//  read_intervals_by_type
// 
// 		read into map<string, t_discrete_intervals>
//________________________________________________________________________________________
void read_intervals_by_type(t_discrete_intervals_by_type&		intervals_by_contig_by_type, 
							std::istream&						f_intervals,
							bool								inclusive_coordinates)
{
	// adjust values to convert from inclusive to open close coordinates
	unsigned cnt_line = 0;
	unsigned adjust_value = inclusive_coordinates ? 1 : 0;

	sp_genomic_interval_by_type_grammar genomic_interval_by_type_grammar = create_genomic_interval_by_type_grammar(cnt_line);

	string str;
	while (std::getline(f_intervals, str))
	{
		cnt_line ++;
		if (!str.length() || str[0] == '#')
			continue;

		t_genomic_interval_by_type interval;
		if (!parse_genomic_interval_by_type(genomic_interval_by_type_grammar, str, interval))
			break; 

		// convert to open / closed by adding 1
		intervals_by_contig_by_type[interval.interval_type].intervals_by_contig[interval.contig].push_back(t_interval(interval.beg - adjust_value, interval.end));
	}

	//	
	//  initialise all
	// 
	for (i_discrete_intervals_by_type ii = intervals_by_contig_by_type.begin(), end = intervals_by_contig_by_type.end(); ii != end; ++ii)
		ii->second.init();
}


//________________________________________________________________________________________
//
//  write_intervals
// 
// 		read into vector<string, t_discrete_intervals>
//________________________________________________________________________________________
void write_intervals(	std::ostream& 					f_intervals, 
						const t_discrete_intervals& 	intervals, 
						bool 							inclusive_coordinates,
						const std::string&				interval_type_name
						)
{
	// adjust values to convert from inclusive to open close coordinates
	unsigned adjust_value = inclusive_coordinates ? 1 : 0;

	// get sorted contig names
	vector<string> sorted_contigs = intervals.get_sorted_contigs();

	namespace karma = boost::spirit::karma;
	using boost::spirit::karma::uint_;
	using boost::spirit::karma::lit;
    using boost::spirit::karma::generate;

	std::ostreambuf_iterator<char> sink(f_intervals);

	BOOST_FOREACH(const string& contig_name, sorted_contigs) 
	{
		BOOST_FOREACH(const t_interval& interval, intervals.intervals_by_contig.find(contig_name)->second) 
		{
			generate(sink,
						(
							karma::string 
							<< karma::string 
							<< lit('\t') << uint_
							<< lit('\t') << uint_
							<< lit('\n')
						),
						interval_type_name, contig_name, interval.beg + adjust_value, interval.end
					);
		}
	}
}

//________________________________________________________________________________________
//
//  write_intervals_by_type
// 
// 		read into vector<string, t_discrete_intervals>
//________________________________________________________________________________________
void write_intervals_by_type(	std::ostream& 					f_intervals, 
								std::map<string, unsigned>& 	map_interval_type_to_index, 
								std::vector<string> 			interval_type_names, 
								t_vec_discrete_intervals& 		intervals_by_contig_by_type,
								bool 							inclusive_coordinates)
{
	using boost::range::sort;
	using boost::range::copy;

	sort(interval_type_names);
	BOOST_FOREACH(const string& interval_type_name, interval_type_names) 
	{
		unsigned interval_id = map_interval_type_to_index.find(interval_type_name)->second;
		write_intervals(	f_intervals, 
							intervals_by_contig_by_type[interval_id], 
							inclusive_coordinates,
							interval_type_name + "\t");

	}
}


//________________________________________________________________________________________
//
//  read_intervals_by_type
// 
// 		read into vector<string, t_discrete_intervals>
//________________________________________________________________________________________
void read_intervals_by_type(std::istream& 						f_intervals, 
							t_vec_discrete_intervals& 			intervals_by_contig_by_type, 
							std::map<string, unsigned>& 		map_interval_type_to_index, 
							std::vector<string>& 				additional_interval_types, 
							bool 								inclusive_coordinates)
{
	// adjust values to convert from inclusive to open close coordinates
	unsigned adjust_value = inclusive_coordinates ? 1 : 0;
	unsigned cnt_line = 0;

	sp_genomic_interval_by_type_grammar genomic_interval_by_type_grammar = create_genomic_interval_by_type_grammar(cnt_line);

	// assign as many copies as necessary for each interval_type in map
	intervals_by_contig_by_type.resize(map_interval_type_to_index.size());

	string str;
	while (std::getline(f_intervals, str))
	{
		cnt_line ++;
		if (!str.length() || str[0] == '#')
			continue;

		t_genomic_interval_by_type interval;
		if (!parse_genomic_interval_by_type(genomic_interval_by_type_grammar, str, interval))
			break; 

		unsigned interval_type_index = 0;
		// additional type
		if (!map_interval_type_to_index.count(interval.interval_type))
		{
			map_interval_type_to_index[interval.interval_type] = intervals_by_contig_by_type.size();
			intervals_by_contig_by_type.push_back(t_discrete_intervals());
			additional_interval_types.push_back(interval.interval_type);
		}
		interval_type_index = map_interval_type_to_index[interval.interval_type];

		// convert to open / closed by adding 1
		intervals_by_contig_by_type[interval_type_index].intervals_by_contig[interval.contig].push_back(t_interval(interval.beg - adjust_value, interval.end));
	}

	//	
	//  initialise all
	// 
	using boost::range::sort;
	for (i_vec_discrete_intervals ii = intervals_by_contig_by_type.begin(), end = intervals_by_contig_by_type.end(); ii != end; ++ii)
		ii->init();
}

unsigned t_discrete_intervals::total_count() const
{
	unsigned cnt_intervals = 0;
	BOOST_FOREACH(const t_intervals_by_contig::value_type& p, intervals_by_contig) 
		cnt_intervals += p.second.size();
	return cnt_intervals;
}

unsigned t_discrete_intervals::total_size() const
{
	unsigned total_interval_size = 0;
	BOOST_FOREACH(const t_intervals_by_contig::value_type& p, intervals_by_contig) 
		BOOST_FOREACH(const t_interval& interval, p.second) 
			total_interval_size += interval.end - interval.beg;
	return total_interval_size;
}


std::vector<std::string> t_discrete_intervals::get_sorted_contigs() const
{
	using boost::range::sort;
	using boost::range::copy;
	vector<string> sorted_contigs;
	copy(intervals_by_contig | map_keys, back_inserter(sorted_contigs));
	sort(sorted_contigs, t_sort_by_embedded_numbers());
	return sorted_contigs;
}

}

